In this report, we examine the impact of a scabies outbreak on antiparasitic treatment trends and the preferred treatment options across Scotland during the second half of 2023 (June to December). The analysis focuses on prescription data to identify patterns and shifts in medication choices among health boards during this period.

Loading Required Packages and Data

First, we will load the different R packages required for this report. These packages are necessary to clean, manipulate, visualize, and analyze the data. Next, we will load the data for the months we are interested in. In this case, we have used a URL to fetch the data, ensuring that the process can be replicated by anyone running the code.

pacman::p_load(tidyverse, janitor, gt, gtExtras, purrr, ggrepel, sf, patchwork, here, leaflet, htmlwidgets, htmltools)
options(timeout = 1000)
urls <- c("https://www.opendata.nhs.scot/dataset/84393984-14e9-4b0d-a797-b288db64d088/resource/66e9611a-fbea-45b2-9edd-351c388fd06d/download/pitc202306.csv","https://www.opendata.nhs.scot/dataset/84393984-14e9-4b0d-a797-b288db64d088/resource/7bb45ee6-6f1c-45f4-933b-958bdbe0ca4f/download/pitc202307.csv","https://www.opendata.nhs.scot/dataset/84393984-14e9-4b0d-a797-b288db64d088/resource/72ca4ad2-0228-4672-9eb0-cc911e4a8ca7/download/pitc202308.csv","https://www.opendata.nhs.scot/dataset/84393984-14e9-4b0d-a797-b288db64d088/resource/2d3d240f-1467-4a91-9f0e-769745650cb9/download/pitc202309.csv","https://www.opendata.nhs.scot/dataset/84393984-14e9-4b0d-a797-b288db64d088/resource/94134438-db1d-4b8d-8f70-5e9b0e47bd03/download/pitc202310.csv","https://www.opendata.nhs.scot/dataset/84393984-14e9-4b0d-a797-b288db64d088/resource/21d11fed-a494-4e30-9bb6-46143c3b9530/download/pitc202311.csv","https://www.opendata.nhs.scot/dataset/84393984-14e9-4b0d-a797-b288db64d088/resource/00cdf8d7-f784-4e87-894c-34f2540ea6ab/download/pitc202312.csv")
data_list <- map(urls, ~ read.csv(.x) %>% clean_names())
dataset_year23 <- setNames(data_list, c("June", "July", "August", "September", "October", "November", "December"))

Loading Additional Data

Next, we will load additional data necessary for the analysis. In this case, we will load a list of the Health Boards and their corresponding names from which we will extract only the Health Board code and the Health Board name. Next, we will load the population data, which contains the population numbers for each Health Board. Since this is aggregated data, we will filter the dataset to select only the rows for “All”. Finally, we will create a list of all the datasets for further processing.

# Loading NHS Health Boards Data
nhs_board <- read.csv(here("data", "Health_Boards_(Dec_2020)_Names_and_Codes_in_Scotland.csv")) %>% clean_names()%>%select(hb20cd, hb20nm)%>%rename(hb = hb20cd, hb_name = hb20nm)
# Loading Population Data
population_data2023 <- read.csv(here("data","hb2019_pop_est_14102024 (1).csv")) %>%clean_names() %>% filter(year == "2023" & sex == "All") %>% select(hb, all_ages)
#Data for GP Practice Codes
gp_scot <- read.csv(here("data", "/practice_contactdetails_oct2023-open-data (1).csv"))%>% clean_names() %>% select(practice_code, data_zone, hb)
#Spatial data with of data zones
data_zones <- suppressMessages(st_read(here("data", "SG_DataZone_Bdry_2011.shp")))
## Reading layer `SG_DataZone_Bdry_2011' from data source 
##   `C:\Users\beaba\Desktop\Neuroscience\Year 4 UoE\Data Science\ICA\B203192\data\SG_DataZone_Bdry_2011.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 6976 features and 9 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 5513 ymin: 530252.8 xmax: 470323 ymax: 1220302
## Projected CRS: OSGB36 / British National Grid
DataZone_with_hb <-read_csv(here("data","DataZone2011lookup_2024-07-04 (1).csv")) %>%clean_names() %>%  select(dz2011_code, hb_code)

Data Analysis

We created a factor including the different treatments for scabies: topical (permethrin cream, malathion lotion, benzyl benzoate lotion, crotamiton cream and sulfur ointment), oral (ivermectin).

Trend Analysis over 2023

To visualize the trend of the scabies outbreak during the second half of 2023 (June to December), we combined all the monthly datasets into one using a list. This allowed us to create a comprehensive dataset that captured scabies treatment prescription trends over the specified months.

By joining the datasets, we ensured that the data was correctly aligned by month and treatment, and could be further analyzed to identify trends, patterns, and differences across health boards.

#Factor with the treatments
scabies <- c("PERMETHRIN 5|PERMETHRIN_CR|MALATHI|BENZYL BENZOA|IVERMECTIN 3")
#Creating a full data set with all the months.
total_data23_sc <- dataset_year23 %>%
  imap_dfr(~ .x %>%filter(str_detect(bnf_item_description, scabies) & !is.na(bnf_item_description)) %>% group_by(hbt, gp_practice) %>%summarise(paid_quantity = sum(paid_quantity, na.rm = TRUE), .groups = "drop") %>%mutate(month = .y))
with_names_sc <- total_data23_sc %>% full_join(nhs_board, join_by(hbt == hb)) %>% mutate(month =factor(month, levels = c("June", "July", "August", "September", "October", "November", "December"),ordered = TRUE)) %>%arrange(hb_name, month)
#Trend over the months 
no_gp <- with_names_sc %>% filter(!is.na(month)) %>% group_by(month, hb_name, hbt) %>% summarise(paid_quantity = sum(paid_quantity))
rate_sc <- no_gp %>% full_join(population_data2023, join_by(hbt == hb)) %>% mutate(rate_per1k = paid_quantity/all_ages * 1000)

sc_squares <- rate_sc %>%filter(!is.na(month)) %>% ggplot(aes(x = month, y = rate_per1k, colour = hb_name))+  geom_line(aes(group = hb_name), size = 1) + geom_text_repel(data = . %>% arrange(desc(month)) %>% group_by(hb_name) %>% slice(1),aes(label = round(rate_per1k, 1)),nudge_x = 0.5,  size = 3, direction = "y",  box.padding = 0.3,point.padding = 0.3, show.legend = FALSE) +geom_point(data = data.frame(hb_name = unique(rate_sc$hb_name)), aes(x = NA_real_, y = NA_real_, colour = hb_name), shape = 15, size = 5,show.legend = TRUE) +  labs(title = "Scabies Outbreak in Scotland (June to December 2023)", subtitle = "Monthly Trends in Scabies Treatment Prescriptions",   x = "Month", y = "Quatity of Antiparasitic Prescriptions /10,000 people", colour = "Health Board") +theme_bw() + guides( colour =guide_legend(title.position = "top", title.hjust = 0.5, override.aes = list(shape = 15,     size = 5)))+theme(legend.position = "bottom", legend.direction = "horizontal") + scale_color_manual(values = c("#e6194b", "#ffe119", "#B0C4DE", "#f58231","#911eb4", "#46f0f0","#f032e6", "#bcf60c","#3cb44b", "#fabebe", "#0000FF", "#e6beff","#9a6324", "#606060"))

sc_squares######check it because the last 2 nhbs cant be seen well

Analysis of Preferred Scabies Medications

In this section, we focus on identifying the preferred medications for scabies treatment across Scotland’s health boards. To do so, we first isolated the relevant data, including the health board, general practice, medication description, and the number of prescriptions for each treatment. Next, we merged the dataset with the list of health boards to label each observation accordingly. This allowed us to associate each prescription with the We will now analyse what are the preferred medication for scabies treatment with a bar plot.

separated_data <- dataset_year23 %>%imap_dfr(~ .x %>%select(hbt, gp_practice, bnf_item_description, paid_quantity) %>%filter(str_detect(bnf_item_description, scabies)) %>% mutate(month = .y))
sep_with_names_sc <- separated_data %>% full_join(nhs_board, join_by(hbt == hb)) %>%filter(!is.na(paid_quantity))
sep_with_population_sc <- sep_with_names_sc %>%full_join(population_data2023, join_by(hbt == hb)) 
sep_bar <- sep_with_names_sc %>% ggplot(aes(x = hb_name, y = paid_quantity))+ scale_fill_brewer(palette = "Spectral")+geom_col(aes(fill = bnf_item_description), width = .5)+theme_classic()+theme(axis.text.x = element_text(angle = 65, hjust = 1))+labs(title = "Preferred Scabies Treatments Across Scotland", subtitle = "Breakdown by health board and treatment type (June - December 2023)", x = "Health Board", y = "Number of Prescriptions", fill = "Antiparasitic Medication")
sep_bar

This bar plot illustrates that the most commonly prescribed medication for scabies was permethrin cream. However, it is important to note that over time, resistance to permethrin has been observed in scabies populations, which may affect the efficacy of the treatment, and thus more affected areas may choose alternative treatments. Additionally, the preference for different scabies treatments can vary between countries, depending on factors such as medication availability and healthcare access.

Relationship Between Different Treatments Over Time

To visualize the relationship between the different scabies treatments and their distribution across months, we have created a table showing the percentage of prescriptions for each medication per health board. We will focus on the health boards that experienced the highest impact, which could depend on the rate of scabies cases in those regions.

The table below presents the proportion of prescriptions for different scabies treatments (permethrin cream, malathion lotion, ivermectin tablets, and benzyl benzoate) from June to December 2023. It also shows how these proportions varied across different health boards in Scotland.

By calculating the percentage of each treatment prescribed relative to the total prescriptions in each health board and month, we can examine trends over time. This analysis highlights the most commonly used treatments for scabies and offers insights into regional variations in medication usage.

This table provides a clear view of the shifts in treatment preferences across health boards, enabling us to identify trends and the relative importance of different scabies treatments. The use of visual tools like this allows for a better understanding of the healthcare response to the scabies outbreak over the selected period.

total_presc_table <- sep_with_names_sc %>%group_by(hb_name, month) %>%summarise(total_prescriptions =sum(paid_quantity)) %>% ungroup()

presc_table <- sep_with_names_sc %>% filter(!is.na(bnf_item_description)) %>%group_by(hb_name, month, bnf_item_description)%>% summarise(treatment_prescriptions = sum(paid_quantity, na.rm = TRUE)) %>% full_join(total_presc_table, by = c("hb_name", "month"))%>% mutate(percentage =(treatment_prescriptions / total_prescriptions) * 100)%>%group_by(hb_name) %>% select(hb_name, month, bnf_item_description, percentage)%>% mutate(month = factor(month, levels = c("June", "July", "August", "September", "October", "November", "December"))) %>%filter(month %in% c("June", "September", "December"), hb_name %in% c("Shetland", "Tayside", "Lothian", "Grampian", "Fife"))%>% arrange(hb_name, month)%>%pivot_wider(names_from = bnf_item_description, values_from = percentage) %>%replace_na(list("MALATHION 0.5% AQUEOUS LIQUID" = 0, "PERMETHRIN 5% CREAM" = 0, "IVERMECTIN 3MG TABLETS" = 0, "BENZYL BENZOATE 25% APPLICATION" = 0))%>% gt()%>%tab_header(title = "Scabies Treatments by Health Board and Month", subtitle = "Data from Scottish Health Boards, June to December 2023") %>% cols_label(hb_name = "Health Board", month = "Month", "MALATHION 0.5% AQUEOUS LIQUID" = "Malathion (%)", "PERMETHRIN 5% CREAM" = "Permethrin (%)", "IVERMECTIN 3MG TABLETS" = "Ivermectin (%)", "BENZYL BENZOATE 25% APPLICATION" = "Benzyl Benzoate (%)") %>% fmt_number(columns = c("MALATHION 0.5% AQUEOUS LIQUID", "PERMETHRIN 5% CREAM", "IVERMECTIN 3MG TABLETS","BENZYL BENZOATE 25% APPLICATION"), decimals = 1) %>% cols_align(align = "center", columns = c("MALATHION 0.5% AQUEOUS LIQUID", "PERMETHRIN 5% CREAM", "IVERMECTIN 3MG TABLETS","BENZYL BENZOATE 25% APPLICATION")) %>% summary_rows(columns = c("MALATHION 0.5% AQUEOUS LIQUID", "PERMETHRIN 5% CREAM", "IVERMECTIN 3MG TABLETS","BENZYL BENZOATE 25% APPLICATION"),fns = list("Average" = ~mean(., na.rm = TRUE)),fmt = list(~ fmt_number(., decimals = 2))) %>% grand_summary_rows(columns = c("MALATHION 0.5% AQUEOUS LIQUID", "PERMETHRIN 5% CREAM", "IVERMECTIN 3MG TABLETS","BENZYL BENZOATE 25% APPLICATION"),fns = list("Overall Average" = ~mean(., na.rm = TRUE)),fmt = list(~ fmt_number(., decimals = 1))) %>%cols_move(columns = "IVERMECTIN 3MG TABLETS", after = "BENZYL BENZOATE 25% APPLICATION") %>%tab_spanner(label = "Topic Treatment", columns = c("MALATHION 0.5% AQUEOUS LIQUID", "PERMETHRIN 5% CREAM", "BENZYL BENZOATE 25% APPLICATION")) %>% tab_spanner(label = "Oral Treatment", columns = c("IVERMECTIN 3MG TABLETS")) %>% tab_style(style = cell_fill(color = alpha("#e6beff", 0.2)),locations = cells_body(rows = hb_name == "Shetland")) %>%tab_style(style = cell_fill(color = alpha("maroon", 0.2)),locations = cells_body(rows = hb_name == "Tayside"))%>%  tab_style(style = cell_fill(color = alpha("lightblue", 0.2)),locations = cells_body(rows = hb_name == "Lothian")) %>%tab_style(style = cell_fill(color = alpha("#D3D3D9", 0.2)),locations = cells_body(rows = hb_name == "Grampian")) %>%tab_style(style = cell_fill(color = alpha("#B4E0B8", 0.2)),locations = cells_body(rows = hb_name == "Fife")) %>%tab_style(style = cell_fill(color = alpha("#FFF2B2", 0.7)), locations = cells_summary(columns = everything())) %>%tab_style(style = cell_text(weight = "bold"), locations = cells_summary(columns = everything())) %>%tab_style(style = cell_fill(color =alpha("#FFB84D", 0.4)),locations = cells_grand_summary(columns = everything())) %>%tab_style(style = cell_text(weight ="bold"), locations = cells_grand_summary(columns = everything()))

presc_table
Scabies Treatments by Health Board and Month
Data from Scottish Health Boards, June to December 2023
Month
Topic Treatment
Oral Treatment
Malathion (%) Permethrin (%) Benzyl Benzoate (%) Ivermectin (%)
Fife
June 41.0 59.0 0.0 0.0
September 17.8 82.2 0.0 0.0
December 10.3 89.4 0.0 0.2
Average — 23.05 76.88 0.00 0.07
Grampian
June 47.1 52.8 0.0 0.1
September 8.9 90.6 0.3 0.2
December 2.6 97.2 0.0 0.2
Average — 19.53 80.21 0.11 0.14
Lothian
June 41.9 58.1 0.0 0.0
September 12.8 87.1 0.0 0.0
December 3.1 96.5 0.1 0.2
Average — 19.27 80.59 0.04 0.09
Shetland
June 69.4 30.6 0.0 0.0
September 16.1 83.9 0.0 0.0
December 0.0 100.0 0.0 0.0
Average — 28.52 71.48 0.00 0.00
Tayside
June 25.2 74.7 0.0 0.0
September 5.3 94.6 0.0 0.1
December 5.6 94.1 0.0 0.4
Average — 12.06 87.80 0.00 0.14
Overall Average — 20.5 79.4 0.0 0.1

In this table, we can clearly see how permethrin is the preferred treatment for scabies followed by Malathion. This results are interesting because the WHO recommends ivermectin treatment for scabies outbreaks… However, in the UK ivermectin was unlicensed after the recent COVID-19 pandemic due to false claims about its efficacy as a treatment for COVID and, as a result, the consumption of ivermectin for horses… However, during the scabies outbreak due to the increasing numbers of permethrin-resistant scabies cases, ivermectin had to be licensed again in October 2023.

#Analysis of the most affected localizations Limitation about the data and its localization. The data was log-transformed for a better visualization of the data zones affected by the outbre

The dataset containing GP practice codes does not cover all data zones, likely because not every data zone includes a GP practice. Additionally, the complete list of data zones does not specify the corresponding Health Board for each zone. Therefore, to plot a map of a Health Board divided by data zones, a supplementary dataset is required to link each data zone to its respective Health Board. This supplementary dataset was loaded at the beginning of the report. Furthermore, this analysis is based on spatial data from the 2011 data zones, which may have undergone changes since then. However, this limitation has been addressed by incorporating the supplementary dataset which accurately maps each data zone to its current Health Board.

with_datangp <- total_data23_sc %>% full_join(gp_scot, join_by(gp_practice == practice_code))%>% full_join(data_zones, join_by(data_zone == DataZone))%>% full_join(DataZone_with_hb, join_by(data_zone == dz2011_code)) %>% mutate(paid_quantity = replace_na(paid_quantity, 0)) %>% select(-hbt, -hb)
with_datangp_sf <-  st_as_sf(with_datangp)



hb_data <- tibble(
  hb_code = c("S08000026", "S08000030", "S08000024", "S08000020", "S08000029"),
  title = c("NHS Shetland", "NHS Tayside", "NHS Lothian", "NHS Grampian", "NHS Fife"),
  lat = c(60.7931, 56.4852, 55.9444, 57.1497, 56.0656),  
  lng = c(-1.3309, -3.0653, -3.1785, -2.1095, -3.1472))

with_datangp_sf <- st_transform(with_datangp_sf, crs = 4326)  # 4326 is the EPSG code for WGS84

# Function to create individual maps without the legend
create_leaflet_map <- function(hb_code, title, data, lat, lng) {
  # Filter data for the specific health board
  filtered_data <- data %>%
    filter(hb_code == !!hb_code)
  
  # Apply log transformation to the paid_quantity for the color scale
  log_paid_quantity <- log1p(filtered_data$paid_quantity)
  
  # Create the color palette
  color_palette <- colorNumeric(
    palette = "Blues",
    domain = log_paid_quantity,
    na.color = "transparent"
  )
  
  # Create the leaflet map without the legend
  leaflet(filtered_data) %>%
    addProviderTiles("CartoDB.Positron") %>%
    addPolygons(
      fillColor = ~color_palette(log1p(paid_quantity)),
      color = "white",
      weight = 1,
      opacity = 1,
      fillOpacity = 0.8,
      popup = ~paste("<strong>", title, "</strong><br>", "Paid Quantity: ", paid_quantity)
    ) %>%
    setView(lng = lng, lat = lat, zoom = 8)
}

# Generate the maps for each health board
maps <- hb_data %>%
  pmap(~ create_leaflet_map(..1, ..2, with_datangp_sf, ..3, ..4))

# Add the image of the legend (after taking a screenshot of the legend in the browser)
legend_image_url <- "data/Legend.jpg"

#C:/Users/beaba/Desktop/Neuroscience/Year 4 UoE/Data Science/ICA/B203192/

legend_image_html <- tags$img(
  src = legend_image_url, 
  width = "300px",    # Adjust width as necessary
  height = "auto",    # Maintain aspect ratio
  alt = "Legend Image" # Alternative text for the image
)

# Combine the maps with the embedded image of the legend
combined_map <- div(
  style = "display: flex; flex-direction: row; justify-content: space-between; flex-wrap: wrap;",
  lapply(1:length(maps), function(i) {
    div(
      style = "margin-right: 20px; width: 250px;",
      h4(hb_data$title[i]),  # Title above each map
      maps[[i]]               # Corresponding map
    )
  }),
  div(
    style = "flex: 1 1 100%; margin-top: 20px; display: flex; justify-content: center;",
    legend_image_html  # Add the image here
  )
)

# Render the combined map with the screenshot legend
browsable(combined_map)

NHS Shetland

NHS Tayside

NHS Lothian

NHS Grampian

NHS Fife

Legend Image
Legend Image
Legend Image

It is worth noting that there are around 300 rows that could be not included because they are under the “gp practice” code 99999, 99998, 99997, 99996, as these are the codes for dentistries, hospitals, unallocated? , and pharmacies. These cant be located in a map.